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Abstract 

This paper focuses on the temporal discretization of the Langevin dynamics, and on different 
resulting numerical integration schemes. Using a method based on the exponentiation of time 
dependent operators, we carefully derive a numerical scheme for the Langevin dynamics, that we 
found equivalent to the proposal of Ermak and Buckholtz, J.Comput.Phys 35, pl69 (1980), and 
not simply to the stochastic version of the velocity- Verlet algorithm. However, we checked on 
numerical simulations that both algorithms give similar results, and share the same "weak order 
two" accuracy. We then apply the same strategy to derive and test two numerical schemes for the 
dissipative particle dynamics (DPD). The first one of them was found to compare well, in terms of 
speed and accuracy, with the best currently available algorithms. 
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1. INTRODUCTION 



Langevin dynamics is widely used in large-scale simulations of colloids, polymers, and 
any system requiring a thermostat, either as a physical ingredient of the system or as an 
expedient in very long simulations to prevent the divergence of the temperature. This 
ubiquity of stochastic dynamics in the realm of numerical simulations has motivated a lot 
of work to design simple, efficient and accurate numerical schemes mimicking the Langevin 
dynamics. For instance, in a recent paper-, Ricci and Ciccotti devised a discrete integration 
scheme for the Langevin dynamics of a set of classical particles. Curiously enough, the 
resulting algorithm differs in a way (apparently) difficult to reconcile with another well- 
known Langevin algorithm, the so-called Ermak algorithnt^ii^. In particular, the Ermak 
algorithm needs two random numbers per degree of freedom for each time step, whereas the 
Ricci-Ciccotti algorithm, which can also be called stochastic Verlet algorithm, requires just 
one for an equal claimed accuracy: as the computational effort is obviously different in each 
case, either there is redundancy in one case, or an incorrect evaluation of the precision of 
the algorithm in the other case. 

One purpose of this paper is to explain where these differences stem from. To this end, 
we revisit the techniques based on Trotter formulas in the context of stochastic dynamics, 
and show that when commutators are taken into account correctly up to the desired order, 
the Trotter formula yields the Ermak algorithm. 

However, we show also that the "simplified" algorithms using just one random number 
per degree of freedom seem actually as good as the correct Ermak algorithm, due to a 
statistical cancellation of the missing term, what the numerical simulations amply confirm. 

The second purpose of this paper follows naturally from these developments: the tech- 
nique based on a Trotter expansion can be easily adapted to the description of the DPD 
dynamics, a stochastic dynamics a la Langevin which moreover preserves the momentum, 
and hence is extremely appealing for applications requiring the stochasticity either as a tool 
to stabilize the dynamics or as an effective "heat bath" in simulations of non-equilibrium 
stationary states^. As this dynamics was quite recently introduced, not many algorithms 
have been proposed for the DPD^"^, some of them having been proved flawed (see for instance 
the discussion in Ref.-). We devise here a new DPD algorithm quite simple to implement 
and which has an efficiency comparable to the Shardlow algorithm^. 



2 



The paper is organized as follows: in the next section, we point out the discrepancies 
which exist between the stochastic Verlet and the Ermak versions of Langevin algorithms. 
Then, in Section [31 we explain how to conduct a Trotter expansion in our Langevin context. 
Concluding this part, we show why two random numbers per step and per degree of freedom 
are, in principle, required in a simulation. The Langevin dynamics algorithms resulting from 
our Trotter expansion are given in Section HI In a second time, we apply our developments 
to the generation of DPD algorithms (Section [5]) and discuss their efficiency in Section [61 
We show in particular that the "bare" efficiency of a DPD algorithm must be carefully 
disconnected from its ability to handle the intrinsic singularity of the dissipation term when 
two particles cross each other (a situation likely to happen only in a perfect gas or with 
particles with soft-core potentials). In Section [21 we also show numerics which endorse as 
safe the use of stochastic velocity- Verlet algorithms, such as the one proposed by Ricci and 
Ciccotti, in a Langevin simulation. 



2. THE TWO LANGEVIN DYNAMICS DISCRETIZATION SCHEMES 

The Langevin dynamics of a single particle with phase space coordinates (position and 
momentum) r{t) and p{t) reads: 

r = - ; P = /(r(t)) - 7P(t) + o-V(t), (1) 
m 

where / is the external force acting on the particle (neighboring particles and/or fields), 7 
a friction coefficient (homogeneous to a frequency), T the temperature of the thermostat, 
cr = y/2m'ykBT, and V'(^) ^ white noise with unit variance (equivalently written as the 
derivative of a Wiener process) : {tpait)) = 0, {ipa(t)ip/3{t')) = 5a^p5{t — t') (a and [3 refer to 
the space indices). 

In one version of Ermak scheme^, velocities and positions are updated according to the 
rule: 



1 - e-^^* p{t) At2 
7 m 2m ' 



r{t + At) = r{t) + ^ + — /(r(t)) + Ar^; (2) 



Pit + At) = e-^'^'p{t) 



/(r(t)) + /(r(t + At)) 



+ Ap^, (3) 



27 

where p(t) and r(t) are the particle's impulsion and position at time t. In the absence 
of stochastic terms Ap^ and Ar^ , and for vanishing friction coefficient 7, this algorithm 



reduces to the well-known velocity- Verlet algorithm. 

By integrating directly ([1]) between and At, and neglecting terms of order higher than 
At^, we obtain the correlations between the two Gaussian increments Ar^ and Ap-^ : 



(Ar^) 
<Ap^> 
ArfArf) 

(ArfApf) 
(ApfApf) 






kT 

kT 

7 

mkT 



27 At - 3 + 4e-^^* - e'^^^' 
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(4) 



1 - e-^^* 
1 - e^^^* 



6. 



a,l3 



(5) 



It is worth noting that Ar^ is of order At^/^ whereas Ap^ is of order At^/^. 
As regards the numerical scheme obtained by Ricci and Ciccotti^, 

At At^ 
r(t + At) = r(t) + — e^^^*/VW + — /(r(t)) + Ar«; 
m 2m 

Pit + At) = e-^^'p{t) + ^/(r(t + At))e-^^*/^ + ^/(r(t))e-3^^*/^ + Ap«, (6) 

it is equivalent to ([2]) and ([3]) up to order At^, except for the random increments Ar^ and 
Ap^ which show slightly different correlations : 



(Ap^) = 
(Ap^ApJ) = a'e--^'At6^,, 
Ar^ - 



(7) 



2m 2m 



In ([7]), the increment Ar^ is proportional to Ap^, and the algorithm requires only the 
generation of a single random number per step. 



Because both algorithms aim at reaching the same accuracy, or at least the same order in 
a power expansion in Af^, it is important to understand the difference between the statistics 
of the random increments. To this end, a natural tool is the Trotter splitting technique, as 
it was invoked in Ref.-. 
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3. MATHEMATICAL PREAMBLE 



3.1. Algebraic formulation of the Langevin dynamics 

The road to finding efficient Brownian numerical algorithms, passes by reformulating 
this dynamical problem in the language of operators. Once expressed in this language, 
the mathematics provide a number of convenient rules, such as the Trotter, or the Baker- 
Campbell-Hausdorff (BCH) formulas, to handle these non-commuting operators. 

This process can be illustrated on a single particle autonomous equation of motion, with 
a mass m set arbitrarily to 1, without loss of generality: 



The resulting dynamics can be subject to a double point of view. The most usual one 
consists in following the phase-space trajectory {r(t),p{t)) of the particle as a function of 
the time t and the initial condition {rQ,po). Drawing an analogy with fluid mechanics, this 
corresponds to a kind of Lagrangian point of view, where the dynamics consists in a mapping 
of the phase-space on itself: 



commonly known as the "Hamiltonian flow". The mapping $ is a time evolution operator 
acting on the phase-space coordinates. This mapping is sometimes loosely symbolized as 
(r(t), = e-^*(r(0), p(0)) but this notation can be misleading as L acts on the phase space 
coordinates, and must not be confused with the Liouvillian operator introduced below. 

The other, dual, point of view consists in focusing on the observables p(r, p; t) of the 
phase space, for a fixed reference point {r,p). This approach is essentially similar to the 
Eulerian description of a flow. For instance, if we consider a probability density function 
(pdf) p{r,p; t) describing the motion of a particle subject to the dynamics ([H]), then the pdf 
taken at two different times t = and t are connected by the relation : 




(8) 




(9) 



p{r,p;t) = p(r',p';0). 



(10) 



where {r',p') is the reciprocal point of {r,p) by the mapping $: 




(11) 



5 



Moreover, for a stationary dynamics, the inversion of $ is just a time reversal, leading 
to {r',p') = $(r,p;— t). Specifying ffTOl) to an infinitesimal time step dt, we obtain the 
Liouvillian dynamics of the observable p{r,p;t) 

with C = p-d/dr + f -d/dp, the Liouvillian operator. The partial differential equation f[T^ 
can be integrated with respect to time, and the pdf p(r, p; t) expressed in terms of its initial 
condition: 

p(r, p- 1) = exp{-tC)p{r, p, 0). (13) 

The time evolution operator exp(— 1£) is a well defined mathematical object, belonging to an 
operator algebra acting on the Hilbert space of differentiable functions of phase space, such 
as p{r,p] t). Basically, exp{—tC) carries the same amount of information as the Hamiltonian 
flow, and its knowledge amounts to the knowledge of the mapping $. An intuitively appeal- 
ing choice for p is the narrow, bell shaped function p^ = exp(— [(r — r^Y + (p — j>o)^]/2e^) 
of width e. Plotting p^{^{r,p; —t)) as a function of r and p for various successive times will 
show the time evolution of a "packet" of independent particles with initial position close 
to (ro,po); each one following its own Hamiltonian trajectory. In higher dimensional phase- 
spaces with deterministic chaos, the packet spreads progressively and eventually occupies all 
the permitted subspace of positions compatible with conserved quantities (ergodicity). The 
e — > limit of this corresponds to an initial condition p{r,p, 0) = 6{r — ro)(5(p — po)] it is 
easy to verify that p{r,p, t) = 6{r — r{t))6{p — p(t)) , where r{t) and p(t) are the solutions of 
the equations of motion with initial condition (rcPo); this shows the complete equivalence 
between the "coordinate" and the "observable" version of the dynamics. 

In our derivation of discrete integration schemes for the Langevin dynamics, we adopt 
this second point of view, consisting of operators acting on the functions of phase space, and 
we stress upon the importance of keeping the minus sign when defining the time evolution 
operator in f|T3|) . Provided this is correctly taken into account, many convenient algebraic 
techniques apply, while getting the same results with the point of view of flows acting on 
coordinates would simply be awkward. Moreover, this convention resembles the familiar 
operator calculus of quantum mechanics (the observable p{r,p;t) playing the role of wave- 
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function). We leave to appendix |X] a simple example showing why working on observables 
gives better results than working with coordinates. 



3.2. The Velocity- Verlet algorithm 

For a stationary dynamics, the next step is to split the time interval t into n small steps 
At and write: 

exp(-t£) = exp(-At£)'^, (14) 

which simply reduces the dynamics to the iteration of a large number of identical elemen- 
tary steps, each one representing the dynamics on a time step At. Finding a good inte- 
gration scheme is equivalent to finding a good estimate of exp(— At£) = exp{A + B), with 
A = —At f{r) ■ d/dp and B = —Atp ■ d/dr. Like in many cases, the complexity of the dy- 
namics comes from the fact that for two non-commuting operators A and B, exp(y4 + B) is 
not exp{A) exp{B), and even though exp(yl) and exp{B) are known, the exponential of the 
sum exp (A -|- B) cannot be calculated. To overcome this difficulty, the strategy is to reduce 
the norm of A and B until their commutator [A, B] = AB — BA becomes small enough to 
apply the Baker- Campbell-Haussdorf (BCH) formula, stating that: 

exp(A) exp(i?) = exp(^A + B + ^[A, B] + ^[A, [A, B]] + ^[B, [B, A]] + , (15) 

with a remaining R consisting of a sum of commutators involving more than four operators 
A OT B each. In the present case, both A and B scale as At, which ensures that for a time 
step small enough, R is just a negligible perturbation term. In practice, the BCH is reversed 
to take the form 

exp(A + B)= exp(A) exp(i?) exp A] - ^[B, [B, A]] + ^[A, [A, B]] + R'^ , (16) 

and one checks that 

exp{A + B)= exp(A/2) exp(5) exp(A/2) exp (-^[5, [B, A]] + ^[A, [A, B]] + R!'^ , 

(17) 
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with R' and R" of order At^. Thus, the two following relations hold, up to terms of order 

exp(-At£) = exp ("^P " |:) exp (-At/(r-) ■ ^ exp ("^P" |:) ; (18) 
- exp . I) exp (-Atp • ^ j exp (-f /(.) ■ |) . (19) 

As recalled in appendix [Cl both splitting are equivalent to the Verlet algorithm. 



3.3. The case of a non stationary dynamics 

When the dynamics is not stationary, the Liouvillian is time-dependent, and the evolution 
operator can be expressed as a time-ordered, or chronological exponential: 

- dsC{s)\=l- dsjC{s)+ dsi ds2 jC{si)jC{s2) . . . 
Jti / Jti Jti Jti 

The symbol T ensures that all the time dependent operators are written from right to left 
in chronological order. Chronological exponentials must be used as soon as the differential 
operators jC{s) and jC{s') taken at two different times, are found to be non-commuting. This 
situation indeed arises when considering the Langevin dynamics of equation ([T]): 

r = -; p = f{r)-jp + ai/^{t). (21) 

m 

The equation for the momentum is modified by a friction coefficient — 7P and a random (or 
Langevin) force il^(t), obeying the usual fluctuation dissipation relation. Straightforwardly, 
we associate to this stochastic equation a "Liouvillian" : 

m = -■!- + (fir) -1P + <Jm) ■ ^, (22) 
m or op 

where the time dependence lies in the random force. By doing so, we implicitly assume a 
Stratonovitch interpretation of the stochastic equation, the only one for which the stochastic 
equivalent of eq. (fT2l) holds^. 

The handling of T exp(— J^*"*"^* i2(s)) shows important differences with the au- 
tonomous, time-independent Hamiltonian situation. First, one cannot simply replace 
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Texp(— J^*^^* by an ordinary exponential without possibly introducing an error of 

order At^/^. Second, when splitting the Liouvillian into A = Atp/m - d/dr and B = C — A, 
the term B is of order At^/^, due to the presence of the random force. This entails that 
the triple commutator [B, [A, B]] in the BCH formula can be of order At^. It is essential 
to either check that this commutator vanishes (as for the usual Langevin dynamics), or to 
keep it in the calculation (as for instance in the DPD situation). 

The basic tool for handling the T-exponential is the Magnus expansion^i, which in the 
present case reads: 

/ j-t+At ^ \ ( j-t+At 

exp I / ds t{s) I - Texp I — / ds C{s] 

exp ( - - / ds du[C{u),C{s)] 

-, pt+At pS pS \ 

--J^ ds duj^ dv[C{v),[C{u)X{s)]]y (23) 

We provide in appendix [B] a simple demonstration of (1231) . It is clear that the right hand 
side reduces to the operator identity for a time independent Liouvillian. Using ( 123|) . one can 
replace the chronological exponential by a simple exponential, which can itself be simplified 
using a BCH formula with triple commutators (fT6l) . 

Finally, let us notice that rather than dealing directly with the stochastic differential 
equation (s.d.e), it is also possible to go to the Fokker-Planck (FP) formulation of the 
stochastic process. In that case, the dynamics expresses naturally in terms of Fokker-Planck 
operators, for which a similar BCH decomposition strategy can be done^, and a numerical 
scheme can be obtained^. In the FP approach, there is no time- dependent term in the 
evolution operator and it is not necessary to consider the time commutators appearing 
in (!23l) . but any connection with the underlying stochastic noise is lost in the process. 



3.4. Application to Langevin dynamics 

We now focus in detail on ( 1211) . The evolution operator can be rewritten as 

/ r-t+At 

Texp I — / ds £(s) 

/ j-t^At \ / r-t+At \ / r-t+At \ 

expi-J^ dsC{s)\ -exp ij^ ds C{s)] ■ T exp i- dsC{s)]. (24) 
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with the product of the two last terms given by the identity (l23l) . Fortunately in this case, 
the triple commutators give rise only to terms of order At^/^, which can be neglected as we 
aim only at a At^ accuracy. Thus, 

exp I — / ds £(s) 1 



"t+At 

Texp ( — / ds C(s\ 



X exp 



ds / du[C{u),C{s)] 

'0 



(25) 



Integrating by part, we get 



t+At 



ds / du[C{u),C{s)] 



a 



At '■*+^* 



t+At 



]_d_ 

m dr 



7 



d_ 

dp 



(26) 



'ip{u)du— / ds / il){u)du 
H Jt Jt 

This time- commutator accounts precisely for the difference between the algorithm of Ricci 

and Ciccotti and the algorithm of Ermak. To support this claim, we simply compare up to 

order At^, the exact time-evolution operator Texp(— Jj*"*^^* ds£(s)), and the corresponding 

simple exponential exp(— j/^"^* ds£(s)), in the case of the s.d.e ([I]). A direct calculation, in 

the exact case, at the order At^, gives: 



r{t + At) = r{t)+p{t) 



At 



m 



f{r{t))-^p{t) 



At^ a '•^^ 



2m m 



ds / il^iu) du\ 



Pit + At) = Pit) + (/(r(t)) - 7P(t)) At + ■ V./(r(t)) - 7/(r(t)) + 7V(t) 



At2 



-a 



t+At I't+At 

il^iu) dw — 7 



f 1 

ds / V'(M)dM>, (27) 
't Jt Jt J 

and is equivalent to applying fl25l) . and to Ermak. The simple exponential case 
exp(— jl^^^ ds£(s)) amounts to replacing the Langevin force by a stepwise averaged quan- 
tity, and to redefining the force accordingly: 

H+At 



fir) ^ fir) = fir) 



a 
At 



xpis)ds. 



It gives: 



rit + At) = rit)+pit)-+{firit))-,pit))— + -- . 



ipiu) du; 



Pit + At) = Pit) + (/(r(t)) - 7P(t)) At + ■ V./(r(t)) - 7/(r(t)) + ^^t) 



(28) 



(29) 

At2 



ct+At pt+At 

+cr / ?/?(«) dw — 70" — / il^iu)du. 



(30) 
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Expression (130|) is equivalent to the stochastic Verlet or Ricci and Ciccotti result. Expres- 
sions fl27|) and fl30l) differ by a term 0(At'^/^), which turns out to be equal to the quantity 
appearing in 



/.t+At /-t+At /.s -1 1-1+ At / 

— y i/>(m) dw — y ds y jp{u)du = J is '0(s)ds. 



(31) 



We note that, in order to get the correct result by using the usual commutator algebra, it 
is absolutely necessary to use the correct convention for the Liouvillian, i.e. exp(— 1£) for 
the forward dynamics of the observables. The above equality is obtained by integrating by 
part the Wiener process. Thus, the Gaussian increments arising in 0271) can be split in two 
statistically independent parts: 

A , 1-1+ At nt+At / A , \ 

^'"-™t/ '^w^'^+^i (t-O^W''^^ 

Ap'' = ff (l - 2^) I'*''' i,{s)ds - laj'*''' (^ - '.') *M<ls, (33) 
each part having zero mean and correlations: 

-t+At \2\ 

^l^{s)ds] j = At; (34) 

t+At 





— IV(s)ds) ) = (35) 



t+At pt+At 



i/,{s)ds I ( ^ - H ip{s)ds ) = 0. (36) 



The statistical independence f l36|) is crucial. It ensures that the difference between (|2j 
and ( |30l) . or equivalently the correction term brought by (!26l) contributes only at the order 
At^ in variance. The current algorithm being a weak order-two algorithm, this term is not 
relevant and negligible at this order. We conclude from this calculation that, when taking 
properly into account all the terms potentially of order At^^^, the algorithm obtained is 
Ermak. However, fl3Tl) turns out to be uncorrelated with the dominant contribution of the 
random force J i/){s)ds, which in turn guarantees that the correction is no longer relevant for 
a weak order-two algorithm. Nonetheless, this statement is not obvious, and the conclusion 
can only be drawn once the corrective terms have been evaluated. 

We conclude by giving a pictorial representation of the difference of these two algorithms. 
From a given noise realization, the stochastic Verlet algorithm gets rid of half of the random 
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numbers that the Ermak algorithm uses, and in a sense, throws away an equal amount of 
information. On the other hand, the generation of the = j {s — At/2)if)[s)ds may seem 
futile, as this is not real information which is "created", but just pure randomness. We 
believe that the issue of this paradoxical situation is that the generation of Ag indeed makes 
sense, even though it is not proper information. To support this claim, we come back to the 
Fokker-Planck (FP) representation of the stochastic process. In the FP representation, one 
aims at calculating the propagator associated to the elementary time interval: 

V{r,p;ro,po;At) = exp^AtCpp) S{p - Po)S{r - tq). (37) 

This function is a bell shaped curve, which is close to a Gaussian distribution. In the two 
dimensional space (r, p), this function is centered around the position that would occupy 
the particle in the absence of noise, and has a "rice grain" shape, with a large axis oriented 
along the p axis, parallel to the vector {At/ {2m), 1 — '~fAt/2). This large axis corresponds 
to the contribution of crAW^ = a j^^^^il^{u) dw. The narrow axis corresponds to the (jAg 
contribution, giving the increment {—aAg/m,a'~fAg) (Fig. [T|, left). By comparison, the 
Ricci and Ciccotti algorithm, in this representation, reduces to a thin line with no width 
(Fig.lll right). No matter how the algorithm is designed, it requires two independent random 
numbers to sample the probability distribution function corresponding to the Fokker-Planck 
propagator. In the same way, when dealing with 3A^ degrees of freedom, the sampling of 
the conditional probability associated to the propagator requires 6A^ random numbers. 

The relevance of the random numbers A^f depends on whether one is interested in a strong 
or a weak convergence of the algorithm. Were the random force ip{t) known, it would make 
sense to estimate the difference between the numerical solution and the exact solution, and 
in that case, the Ermak scheme would be superior to the stochastic Verlet scheme, making it 
a true strong order 1.5 algorithm. However, as far as numerical simulations are concerned, 
we are interested in averaging the observables over many realizations of the noise ip{t) 
(either by running the simulation long enough, or by repeating many times the simulation 
process) and this corresponds to the weak convergence concept, where only the difference 
between the averaged simulated quantities and the averaged exact quantities matters. As 
the variance of A^f is At^, as the variance of AW is only known up to order At^, and as Ag 
and AW are statistically independent, the average contribution of A^f is at best of the same 
order of magnitude as the uncertainty affecting AW. Not providing real "information", the 
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generation of the numbers just makes the Ermak procedure a shghtly better samphng 
of the underlying stochastic dynamics than the other algorithms based on a single random 
number. However, the difference is not quantitative: both schemes are of weak order two. 

As far as the algorithms based on the Fokker-Planck approach are concerned, such as the 
one of Serrano et alM^, discussed in the DPD section below, the resulting algorithms are 
also weak order two. This comes from the fact that the Fokker-Planck approach assumes an 
implicit average of the random noise, when writing the Fokker-Planck evolution operator. 

4. LANGEVIN DYNAMICS OF FLUIDS 

We now specialize the discussion to the Langevin dynamics of a fluid, and consider TV 
particles with usual pairwise interactions. The Liouvillian reads: 

N 3 



^ ' fia{{r}) -'yPia + (Ttl^tait) ^ 



m aria \ J op. 



(38) 



1=1 a=l '- 

where a refers to the space index. As above, the time ordered exponential is replaced by a 
product : 

/ pt+At „ \ / i-t+At ^ \ / ^ pt+At ^ ^ \ 

Texpf — / jC{s)ds j = exp f — / £(s)ds j exp f — - / ds / du [jC{u) , jC{s)] j ; 

= exp(-A-5) -expC, (39) 



with 



I . m fir m fir ' ^ ^ 



1=1 a=l 1=1 a=l 



r-t+At N 3 Q 

B= ds^^[/ia({ri}) -7Pi^ + (T'0^a(s)]^— ; 

i=l a=l '^^^ 
N 3 ^ pt+At ^ 

= HH\^^ [/-({^^}) - wUr)] + a Via(n) du \ 

i=l a=l ^ ^ 



d 



(41) 
(42) 



Using the BCH formula, we obtain two factorizations : 



exp(-A - B) ■ exp(C) = exp - — j ■ exp(-S) ■ exp - — j ■ exp(C); (43) 

/ B\ / B\ 

= exp - -j ■ exp(-A) . exp - -j ■ exp(C). (44) 
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If we forget the commutator term C, we obtain exactly the sphtting proposed by Ricci and 
Ciccotti, but with the opposite sign convention. The commutator term is : 



C 



1 



ds / dn[£(s),i:(n)]; 



[At /■*+^* 

i=l a=l L 



t+At 



ds J il^ia{u)du 



1 d 



1- 



d 



The term exp(C) adds an extra contribution of order At'^/^ to the random Gaussian in- 
crements. For completeness, we give below the correct discretization schemes associated to 
the two splitting ( H3l) and (jBl). A detailed calculation is provided in appendix ([C]). The 
splitting (143|1 leads to 

p^t + At) = pUt) (l - 7At + + fie [r{t) + (l - A^ 



/-t+At pt+At 

+ a i>ia{u) dn - (77 / 

jt Jt Jt 

, . At , , / 7At\ „ , , At^ 
r,^{t + At) = Vi^it) + —pi^{t) ( l-^j 



r 

ds y tl)ia{u)du] 



(46) 



+ 



t+At 



ds / il^ia{u) du. 



(47) 



The random Gaussian increments Ap^, Ar^ appearing in P6|l and P7|) can themselves be 
split in two parts Ap-^ = Ap^ + Ap*^, and Ar^ = Ar^ + Ar'-^. 



cr 



Ar.t = - 



t+At 



ds i>ia{u) du] 



Ar^ 



Apf. 



a [At /■*+^* /■*+'^* 

J '4^ia{u)du- / ds / ii)ia{u)du 

»t+At /-t+At 

cr / V'ia(M) dw - 0-7 



ds ^ V'ia(w) dw; 



At /•*+^* 
cr7 ^ / V'ial^i) du 



t+At 



ds / ij)ia{u)du 



<y\l 



-yAt 



t+At 



du. 



(48) 
(49) 
(50) 
(51) 
(52) 
(53) 



Clearly, Ap'-^ and Ar'^ come from the time commutator (H5l) . It is safe to replace the 
complex random increments Ap^ and Ar^ by their simpler counterpart Ap^ and Ar*^ 
without altering the weak order two nature of the algorithm. 
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The splitting (jHj) leads to another expression: 
pUt + At) = p^t) h - 7At + j + — iU{r{t)) + fUr{t + At))] - ^—f,^{r{t)) 

pt+At pt+At PS 

+cr / il)ia{u) du - •ya / ds ij^ia{u)du] (54) 



.X At , , f At\ „ , , At2 
ri„(t + At) = ri^{t) + —pUt) 1 - 7— + /i«(r(t)) 



ra \ 2 J 2m 

H — / ds i ipia{u) du. (55) 
m Jt Jt 

The above algorithm is equivalent to ([2]), ([3]) and (jlj). The random Gaussian increments are 
obviously identical at the order At^ considered here. Equations (H6l) . (H7j) . ( 15^ and ( l55l) are 
the main result of this section. 

We believe that the absence of Ar'" and Ap'" in Ref.- is due to a time-ordered exponential 
hastily replaced by an ordinary exponential. These terms should have been present, and 
there absence can be traced back to the forward operator J' appearing in the Susuki formula 
exp(Ati7), which, in the presence of the random noise, must acquire some kind of time- 
dependence, such as exp(Atj7'(t)). The derivation of the algorithms of Ref.- is incomplete, 
with a misprint in the result ((51), ([6]). 



5. DISSIPATIVE PARTICLE DYNAMICS 

5.1. Derivation of two original DPD algorithms 

We now apply the same strategy to the dissipative particle dynamics (DPD) introduced 
in Ref.— and discussed in Ref.— li^. The phase space dynamics of the DPD is given by a 
set of stochastic equations: 



Pi. 

7 



(56) 



Pi = fi{{r})+ -—w{rij)'^eij-{p,-pj)e,j+ ^ aw{rij)eijipij{t). (57) 

Each particle i exchanges randomly some momentum with a selected set of neighbors j. This 
exchange dynamics depends on an arbitrary short range function w{rij), and acts along the 
separation rj — Vi, where Vij stands for vj—ri, for \ \rij\ \ and the unit vector e^j for Vij/rij. 
To each pair of particles is associated a random noise ipijit)dt, even though only a fraction 
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of the pairs are truly interacting together at a time, due to the finite range of w{r). The 
prefactor a is equal to \/2^kBT in agreement with the fiuctuation dissipation theorem. Note 
that here, the dissipation coefficient 7 is defined in such a way that the product of 7 and a 
velocity is a force. In order to stay consistent with the previous section, we introduce also a 



friction coefficient 7 = 777,7, homogeneous to a frequency, and such as a = y2rnrfkBT . These 
equations can be rewritten with 3A^-dimensional vectors and matrices: r = p = {pi}, 
dp = {d/dp^}, dr = {d/dr,}, 

RUr-t) = <j\Y. - W \ ; (58) 

(um J 

The dissipative term d{r,p) can be expressed as the product of a space dependent 3A^ x 3A^ 
matrix T{r) and the 3A^-vector momentum p. If a,f3 refer to the space indices, then the 
components of the matrix T[r) read: 

r(r)iaj/3 = 7 2-. —^:2^{ria - rka){rip - rkp){6ij - dkj); (60) 
dir,p) = -r(r)-p. (61) 

It can be checked that the correlations of the noise obey 

{i^^,{t)Mt')) = {s^kSji + 5uSjk)S{t-t'y, (62) 

(Rir-t),^) = 0; 

{R{r-t)i^R{r-t')jfs) = 2-fmkBT6{t - t')T{r)i^j(s. (63) 

To the s.d.e. (|7r|) we associate the Liouvillian 

m =(^^-dr + fir) ■ dp + d(r, p) ■ dp^ + Rir; t) ■ dp. (64) 

The main difference, compared with the ordinary Langevin dynamics, is that the prefactor 
of the noise depends explicitely on r. It is known, however, that this multiplicative noise 
is not dependent of the Ito or the Stratonovitch interpretation of stochastic dynamics^^. 
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Regarding the simple exponential part exp(— J C{s)ds), we define, as usual 



A = At—-dr; 

m 

B = At[/(r) + d(r,p)]-9p + 
and the 3A^ vector K{r;t): 



t+At 



R{r; s)ds 



■d. 



pi 



K{r-t) = f{r) 



At 



t+At 



R{r; s)ds 



(65) 
(66) 

(67) 



A scales like At, B scales like v^At, K{r; t) scales like At One possible splitting is 



Texp 



t+At 



C{s)ds ) = exp ( — — ) ■ exp ( — A] ■ exp 



f) 



■exp(^[i?,[A,B]]-i 



t+At 



ds / du[C{u),C{s)] 
t Jt 

t+At PS PS 

ds du dv[C{v),[C{u),C{s)]]], (68) 



and the other is 

^"t+At 

Texp 



C{s)ds ) = exp ( - ^) ■ exp - 5) ■ exp - ^ 



exp --[B,[A,B]]-- 



t+At 



t 

t+At 



ds / du[C{u),C{s)] 



ds du dv[C{v),[C{u),C{s)]]]. (69) 
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The splitting gives the largest contribution, while the commutators give smaller correct- 
ing terms. The first splitting exp(— ^4/2) exp(— i?) exp(— ^4/2) gives after neglecting terms 
smaller than At^ (algorithm "ABA"), 

At 



r(t + At) = r(t) + 



2m 



p{t + At) 



p(t) +p(t + At) 



At 



(70) 



l-Atr(r-(t) + — p(t) 



.(t) + -p(t) 



Pit) 



+At 



At^f , , At , , 



while exp(— fi/2) exp(— A) exp(— i?/2) gives 

At At^ At^ 

r(t + At) = r(t) + -ip(t) - ^r(r(t)) ■ p(t) + i^K{r{t); t); 

At2 



(71) 



(72) 



p(t + At) = pit) 



At 



At 



r(r(t)) + r(r(t + At)) 



li:(r(t);t) + K(r(t + At);t) 



P{t) 
Ar 



r(r(t))^-p(t) 



r(r(t))-K(r(t);t). (73) 
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We now turn to the double and triple commutator corrections. The correction 
exp(— 1/2 jj duds [C{u),jC{s)]) gives terms proportional to jl^^^{s — At/2)tlJij{s)ds, which 



have been shown to be uncorrelated to the 



j^^^^ipij{u)du , of vanishing mean value and 
of variance At'^. These terms are the pendant of the corrections in the Ermak algorithm, 
and can be neglected in practice. 

After some calculations exposed in appendix ([D]), it appears that together the correction 
of the triple commutators is either vanishing for the splitting fl69|) . or equal to 



a^Ae _d_ 



(74) 



for the splitting 



A popular choice consists of using a triangular shape for w. 
w{r) = l — r/vc, if r < re] 



(75) 



w{r) = 0, if r > Vc- 

This choice ensures that w"^ is differentiable everywhere but in r = 0, and the corrective 
term is well defined. 



I], 



dri, 



-2- ^1 



(76) 



This suggests the final form of our DPD algorithm, corresponding to the splitting ( 1721) and 
SD (algorithm "BAB"). 

(77) 



'n(f'\ At Af^ 

r{t + At) = r{t) + P^At-—T{r{t)).p{t)At + —K{r{t)]^ 

At^ 

r(r(t)) + r(r(t + At)) 



At 



p(t + At) = p{t)- — 



p{t) + ^nr{t)f.p{t) 



+ 



At 



2 

a'^Ae 



K{r{ty,t) + K{r{t + Aty,t) 



Af 



T{r{t))-K{r{ty,t) 



Am 

where the last term g^^ ■ r requires a 3A^ x 3A^ operator matrix 



wiTij) — 

3 ar,„ or 



i/3 



(78) 



(79) 



and 



w{rij) 



(80) 



Formula (1701) and (17T]) . on the one hand, and (177|) . (175]) and (1701) on the other hand, are 
the main results of this section. In the following section, we give more details on how to 
implement these formulas in practice. 
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5.2. Comparison with other existing algorithms 

The BAB algorithm resembles the DPD-velocity Verlet (DPD-VV) algorithnt^. Both 
use an estimate [f{r{t)) + f{r{t + At))]/2 for the forces, and a similar expression for the 
position dependent dissipative forces. Positions are updated in an identical way, positions 
and velocities are updated in the same consecutive order. However, a noticeable difference 
concerns the use of the random numbers : in our case, we use a single random number 
\p{u)du, while DPD-VV use a combination of f^_^^^^i/){u) du and j^^^^^'^il){u) du. 
Moreover some terms differ at the order At^. 

Starting from a Fokker-Planck formulation of the stochastic dynamics, the approach by 
Serrano et alM>^ gives directly a weak order two algorithm, and by-passes all the complica- 
tions related to the time- dependent double and triple commutators in the Magnus expansion. 
Although we do not present a demonstration here, it is also possible to derive their algo- 
rithm starting from our splitting fl69l) . One of the steps in the derivation consists in splitting 
further the exponential of the operator B defined in ( l66ll into a product over all pairs of 
particles. 

In an improvement over DPD-VV, Pagonabarraga et ali^ proposed a self-consistent pro- 
cedure for computing the velocity dependent friction term. This makes the so called self- 
consistent DPD-VV (SC-DPD-VV) a reversible algorithm, for which, in the absence of 
potential forces, the stochastic trajectories can be traced back upon momentum reversal, 
according to the detailed balance probabilities. 

We did not find mention of the reversibility properties of the algorithms of Shardlow S2 
and Serrano et ai, but we believe that, much as the Ermak procedure for Langevin dy- 
namics, the algorithm of Serrano et al. stands a good chance to be reversible, because the 
momentum exchange stochastic dynamics is treated without approximation. The algorithm 
S2 of Shardlow, which almost reduces to the former in the absence of potential forces, can 
only be approximately reversible, while retaining the weak order two convergence property. 

As far as our algorithms are concerned, the stochastic evolution of the momenta p{t), 
with fixed positions r{t), corresponds to a multidimensional Ornstein-Uhlenbeck process, 
as illustrated by eq. (]D4p of appendix [Dl Reversibility is lost when approximations such 
as exp(— FAt) ~ 1 — FAt + are introduced. As keeping exponentials of matrices 

is not algorithmically feasible, the final form of the algorithms BAB and ABA are only 
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approximately reversible, up to the expected order of convergence At^. 



6. NUMERICAL TESTS AND DISCUSSION 
6.1. Ermak versus stochastic Verlet algorithms 

We noticed a fundamental difference between the Ermak and Ricci-Ciccotti (or stochastic 
Verlet) algorithms, the former being in principle slightly more accurate in simulating the 
Langevin dynamics. It is interesting to probe numerically whether such an assertion is really 
detectable. To this end, we considered the simple Ornstein-Uhlenbeck process : 

X + X + X = ipit), (81) 

{m^it')) = m-t'). (82) 

which is exactly solvable. In particular, we have 

(x(t)x(O)) exact = e~*/' X (^cos(V3t/2) + sin(V3t/2)^ . (83) 

A way to check the performances of an algorithm is to measure how fast and how accurately 
it reproduces such a correlator. Thus we define the error associated with a numerical result 

as 



1 



Error = -— / dt 



Tav Jo 



(x(t)x(0))cxact - (t)x(O)) simulated 



(84) 



(more exactly the discrete version of this formula). Tav is a predefined time kept constant. 
This error is a function of the random numbers used to generate (x(t)x(0))sijnuiatcd5 and can 
be considered as a function of the simulation time tsim , as the average (...) is constructed 
by accumulation (running average). Were the time step infinitely small, would such a curve 
decrease toward zero; actually, it saturates (with fluctuations) at a higher value due to 
the unavoidable discreteness of the simulations. Constructed with a single simulation, this 
measure of the error as a function of the simulation length is a noisy curve, and it is more 
relevant to consider the curve obtained by averaging several such errors. 

The result of such a procedure is plotted on Fig. [2], for eighty runs of 10^ time steps (full 
and dashed curves — the light dots show the curves for four subsets of twenty curves to 
show the dispersion). 
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The chosen time Tav = 6 contains the essential part of (x(t)x(O)) (see inset); the time step 
is large : At = 0.1, large enough to reveal in that situation the clear superiority of the Ermak 
algorithm over the Ricci-Ciccotti one : for those large time steps, the Ermak algorithm 
mimics a thermal noise which is less colored than Ricci-Ciccotti's, and hence is a more faithful 
realization of Langevin dynamics. Obviously, when reducing the time step, the spectrum of 
the noise is flatter and flatter for both algorithms and the discrepancies disappear : this is 
shown on the curves in large dots, obtained for a time step At = 0.001 and 10^ simulation 
points (the associated physical time is unchanged); they show that a computational effort 
a hundred times greater than before leads to a rather modest improvement as far as the 
Ermak algorithm is concerned. But the Ermak algorithm can be useful for situations where 
large time steps can be used, without compromising a faithful sampling of the underlying 
Hamiltonian dynamics. We think that these results suggest to always prefer a larger time 
step with the Ermak algorithm to a smaller one with the lightest (in term of computational 
cost) Ricci-Ciccotti scheme. 



6.2. Implementation of the ABA algorithm 

Two algorithms of similar predicted accuracy are available. In practice, after having 
tried them both, we found that the algorithm fl70f71l) was by far the simplest and the most 
efficient of the two. Therefore, we devote to this scheme the largest part of our discussion, 
while briefly mentioning the other at the end of this section. We chose to refer to it as the 
"ABA" algorithm, by reference to the splitting from which it originates. 

It reads explicitly: 

"(a) < = Viit) + ^Mt) 

(b) X, = [r(r')p(t)].At - K,{r')At 

= 7AtE„a^, ^i<rP^M, + ^VA^E„a^, ^K^w., - Mr') At 

(c) = [r(r')X(t)], 

(d) p,(t + At) = pi{t) -X, + -^Yi 

2m 

(e) r,(t + At) = + fp^jt + At), 

where we defined are uncorrected Gaussian random variables 
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with zero mean and unit variance. This algorithm can be made even more efficient by 
merging the first and last steps during the simulation (however, the computation of statistical 
quantities must be performed compulsorily between these steps, otherwise a noticeable bias 
is introduced). The computational effort per time step is slightly higher than the Shardlow 
algorithm (eq. (9) of Ref.- hereafter termed SI algorithm according to the terminology used 
in^), as here two computations implying the sums ^^jW^rij) are required, instead of just 
one for SI. This is the price to pay to have a second order algorithm. Actually we measured 
only an increase of ca. 5% of the time simulation with our algorithm with respect to SI (to 
be efficient, the program must maintain a subdivision of the simulation box and a related 
list of occupation, or alternatively a neighbor list which already exists anyway for MD, such 
that a loop J2j '^{^ij) is fast). 

6.3. Benchmarking with existing algorithms 

When a new algorithm is introduced, it is common practice to rate its performance, and to 
compare it with the other existing algorithms. The usual benchmark includes a calculation 
of the pair distribution function g{r) of a noninteracting fluid (ideal gas), a monitoring of the 
thermalization properties, and a presentation of the self-diffusion data. The absence of a non 
trivial but exactly solvable model makes that there is no stringent test available for assessing 
the ability of an integrator to reproduce faithfully the true DPD dynamics. Previous studies 
discovered that most DPD algorithms do not predict correctly the fiat shape of the g{r) 
characterizing the ideal gas. Instead, some unphysical structure of g{r) near the origin is 
observed, betraying a dynamics with unsatisfactory detailed balance properties. 

We present data for both BAB and ABA algorithms, along with the Shardlow SI and 
DPD-VV algorithms (which are a priori weak order one algorithm). The scale and axis 
of our graphs make our data readily comparable with the one presented by Nikunen et 
air- or Serrano et al.—. Fig. [3] displays the g{r) distribution, for the usual weight function 
w{r) = 1 — r/r^, and three increasing values of the time step At = 0.01, A = 0.05, and 
At = 0.1. Our g{r) is very satisfactory for At = 0.01, comparable with DPD-VV but worse 
than Shardlow SI for At = 0.05, and definitely worse than its competitors for At = 0.1, 
a quite large value by the way. We also plot the BAB algorithm deprived of its correcting 
term, for which the case At = 0.01 clearly demonstrates that this correcting term precisely 



22 



contributes very much to the accuracy of the scheme at small At, and supports firmly the 
validity of our calculation. 

Fig. m summarizes the thermalization properties. The left panel of Fig. H] shows the sim- 
ulated kinetic temperature (T) , compared with the value Tthcrmostat = 1 that the thermostat 
attempts to enforce. In the right panel of Fig. HJ we present a semi-logarithmic plot of the 
difference between the actual simulated temperature and its target value | (T) — Tthcrmostat | , 
for the same three time steps At = 0.01, 0.05 and 0.1. We draw from these data the 
same conclusion as for the g{r) : both algorithms ABA and BAB are found inaccurate for 
At = 0.1, they are comparable with DPD-VV for At = 0.05 while standing below Shardlow 
SI, and finally become as satisfactory as the other algorithms for At = 0.01. 

Fig. El proves that our algorithms predicts correctly the shape of the pair distribution in 
the presence of interaction, for the "model B" introduced in Ref.-, where an interparticle pair 
potential equal to aw^(r)/2; a = 25 is present. Fig. [6] summarizes the results obtained for 
the self-diffusion coefficient D, in the presence of interactions between particles ("model B"). 
In the absence of any exact result for this non-trivial situation, we take the self-diffusion 
value D{At = 0.01) obtained for the smallest time step as a reference, and compare it 
with the values obtained for larger time step At = 0.05 and 0.1. It is also possible to do 
the same study for an ideal gas, but we found that all the tested algorithms (DPD-VV, 
Shardlow SI and ABA) lead roughly to the same value. We believe that the data presented 
here are more discriminant, with a combination of caging effect due to the interactions and 
correlated diffusion due to the non-local dissipation term. 

6.4. The singularity of the weight functions w{r) 

The ideal gas has been widely used as a benchmark to test DPD integrators, because of 
the notorious artefacts in g{r) from which so many DPD algorithms suffer. We claim that 
some attention must be paid in that case: for a choice of a weight function with w{r = 0) 7^ 0, 
the DPD dynamics is singular whenever two particles occupy the same position; therefore, 
the ideal gas (the model par excellence where superposition of particles can occur) should 
only be used as a benchmark model for testing DPD algorithms for weight functions such 
that w{r)/r does not diverge when r ^ 0. Otherwise, what is measured is also the ability of 
the algorithm to deal with this singularity, whereas the algorithms are always written under 
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the assumption of local differentiability of the functions defining the model. This point is 
crucial: a high order numerical scheme may enhance the effect of a singularity and perform 
worse than its lowest order competitors. 

Fig. [7] demonstrates how efficient a regular weight function is at removing the artefacts in 
the pair distribution of ideal gases. It shows the g{r) obtained for Wi{r) = {r /r^) x (1 — r/r^), 
while the thermalization properties of wi{r) are already shown in Fig. HI The weight function 
1/71 (r) vanishes for r = 0, precisely when a pair of particles is colliding, and preventing the 
sudden reversal of the friction forces associated to this pair. Clearly, the use of Wi{r) for soft 
core particles generate no artefact in g{r), without reducing the thermalization efficiency of 
the DPD thermostat. The comparison between the different algorithms turns now to the 
advantage of our ABA scheme, with a flatter g{r) than both DPD-VV and Shardlow SI 
(right panel of Fig. [7]) . 

6.5. Discussion 

First, we would like to emphasize some shortcomings of the above tests, focusing only 
on static correlators : from the measurement of g{r), it is not possible to really check how 
faithful to the true DPD dynamics these algorithms extreme illustration of this 

idea, a Ermak-Langevin algorithm could be qualified as a good DPD algorithm, if only 
static quantities are considered. Actually, a genuinely discriminant test of DPD algorithms 
can be made only in the spirit of the above comparison between Ermak and Ricci-Ciccotti 
algorithms, that is using as a benchmark a dynamical correlator. 

A second question concerns the value of the damping coefficient 7, which for the sake of 
the comparison with previously published data, was fixed to a quite large value 7 = 4.5, 
and for which an efficient thermalization is achieved at the expense of a "sluggish" dynamics. 
Usually, second order algorithms are designed to match smoothly the underlying Hamiltonian 
dynamics in the limit 7^0. Thus, for large values of 7, the putative force of our algorithms 
is not apparent, compared with simpler order one schemes. 

A related question arises : provided that a numerical scheme samples faithfully the 
canonical ensemble and conserves well the total momentum, is that of any importance to 
stick to the actual DPD dynamics or not? The answer is certainly yes if one is interested in 
making connections between numerical results and theoretical developments. In this respect. 
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our second order algorithm is useful, as it matches better Hamiltonian and DPD dynamics 
in the limit 7 — > 0, an interesting limit for those who want the DPD only to thermally 
stabilize a long simulation without perturbing its natural proper dynamics^. 

In our view, it is still an open question to understand why the algorithm of Shardlow 
and of Ref.— are so efficient in avoiding the artefacts of g{r) in Fig. [31 We simply attribute 
this to the sequential update of the velocities over all the pairs of particles, rather than 
doing it all at once, as we propose. We do not, however, consider Fig. [3] as a failure 
of our algorithms, but much more as a warning message for the DPD users interested in 
soft core potentials, concerning the inadequacy of the weight function w{r) = 1 — r/rc in 
that case. Simulations with soft core particles should be done with a regular weight such 
as Wi{r) = (r/rc) x (1 — r/rc). Having tried other possible weight functions, such as a 
parabolic shape 1 — {r/rcY, or the quadratic shape {r /rcY ^ (1 ^ '^l'^c\ we found that the 
thermalization and structure properties are respectively close to the results obtained for the 
singular w(r) and regular W\{r) weight functions, emphasizing the importance of the limit 
value lim^^o w{r)/r, as the key feature for reaping the best of our second order scheme. 

7. CONCLUSION 

We demonstrated that a careful use of a Trotter splitting for Langevin dynamics leads 
to a discretization scheme equivalent to the algorithm of Ermak. The main difficulty in the 
procedure comes from the fact that the "Liouvillian" associated to the Langevin dynamics 
is time dependent, and the contribution of the random noise in the stochastic equation, 
brings significant changes compared with the case of an ordinary differential equation. As 
a consequence, it is necessary to take into account a commutator term given by fl23|) . For 
ordinary Langevin dynamics, this corrective term simplifies to (l25l) . This term is responsible 
for the second independent Gaussian increment appearing in the Ermak algorithm. 

When checking how accurately these numerical schemes can reproduce a known time- 
dependent correlation function, we found that keeping both independent random Gaussian 
increments leads to better results for large integration time steps. However, the benefits of 
the second random number gets smaller when the integration time step is reduced, and in 
the small time step limit, where both algorithms share the same accuracy, a single random 
Gaussian increment give satisfactory results. 
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Using the same method, we were able to propose two algorithms for the dissipative 
particle dynamics (DPD). One of them, termed "ABA", was found to be easy to implement 
and efficient, and we recommend it as a reasonable alternative for applications. Though less 
efficient, the other algorithm, termed "BAB", deserves to be mentioned as an original and 
non trivial weak order two algorithm. We observed that for simulations of soft core particles, 
the unphysical structure of the pair correlation function g{r) is largely a consequence of a 
non vanishing weight function w{r) near r ^ 0. This artefact almost disappears when 
the weight function vanishes at the origin, without weakening the efficiency of the DPD 
thermostat. It is in this limit only that high order algorithms give their best. 

The theoretical developments exposed in this work are poised to provide alternative 
derivations of some existing algorithms. Unlike the Fokker-Planck way, our approach offers 
the possibility to discuss both strong and weak convergence, a fact which extends its scope 
beyond the realm of MD simulations. 

The application to DPD revealed a situation of intermediate complexity, standing be- 
tween the usual additive noise Langevin dynamics, and the most general multiplicative 
noise stochastic equations. In particular, the treatment of a Langevin dynamics with hydro- 
dynamical interactions falls into the same category as DPD, and our discussion will certainly 
help to improve the design of numerical integration schemes, also in that case. 

APPENDIX A: THE SIGN CONVENTION 

To convince ourselves of the importance of adopting the right convention when defin- 
ing the Liouvillian, we consider a simple example. We define the flow 1 as the ordinary 
differential equation (o.d.e): 



X = p; p = 0; 



(Al) 



and the flow 2 as 



X = 0; p = —p. 



(A2) 



Integrating these flows is easy and one flnds: 




(A3) 



and 




(A4) 
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Meanwhile, the LiouviUians associated to these flows are 



ti = p-^ and £2 = -pS-- (A5) 
ox op 

The commutator of Ci and £2 can be calculated, and its value is 

[4,£i] = -P^; (A6) 

We now integrate along the flow 1 for a time interval of At and along the flow 2 for a time 
interval of An. We get the exact result 

$2 (Am) o $1 (At) (x, p) = {x + pAt, pe"^") , (A7) 

which, up to the order 2 in At, Au is equivalent to 

$2(An) o $i(At)(a;,p) = (^x + pAt,p{l - Au + Au'^ /2)^ . (A8) 
Integrating the flows in the reverse order gives a slightly different result 

$i(At) o$2(An)(x,p) = (x + pe-^"At,pe-^"), (A9) 

which, up to the order 2 in At, An is equivalent to 

$i(At) o $2(An)(x,p) = (^x + pAt-pAtAu,p{l-Au + Au'^/2)^. (AlO) 
We observe that when commuting the flows, the following formula holds 

$2(An) o $i(At) = $2i(AnAt) o <l>i(At) o <l>2(An), (All) 
with a flow $21 associated to the o.d.e 

x=p; p = 0. (A12) 

and a Liouvillian £21 equal to 

£21 =P^; (A13) 

In other words, £21 is the opposite of the usual commutator, for operators acting on their 
right hand side: 

£2i = -[£2,£i]. (A14) 
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However, when dealing with the Liouvilhans, the following formula holds: 

e-^"^^e-^*^V(a;,p) = e(-^*)(-^")[^^'^ile-^*^^e-^"^V(a;,p), (A15) 
with, in particular, 

e~^^^^ p{x,p) = p{x —pAt,p) 

^AtAuic^Alp^^ _ pAt,pe^") = p{x - pAt - pAtAM,i)e^"). (A16) 

We conclude that the algebraic formula (lAlSp involves the usual commutator of the differ- 
ential operators (as, for instance, in quantum mechanics), while the equivalent formula ( lAlip 
for the flows require a different definition of the commutator. 

Thus, it is improper to write a flow as an exponential 

$(At) = e^*^, (A17) 

as it would imply 

gAf£igA«£2 ^ g-AMAt[£i,£2]gAM£2gAt£i (A18) 

unless the operators appearing in the Liouvillian are meant to act on their left hand side, 
such as If, fortunately, in most cases the notation (]A17p does not have serious 

consequences, it cannot be used as soon as commutators are explicitely required in the 
calculations. 



APPENDIX B: FROM CHRONOLOGICAL TO REGULAR EXPONENTIALS: 
THE MAGNUS EXPANSION 

We are now trying to evaluate 

exp|^^ £(s)ds^ -Texp 1^-^ £(s)dsj (Bl) 
We split At in n subintervals At/n and write: 

/ ^At \ 'i-l / ;-(j+l)At/n \ 

Texp\^-J £(s)dsj = T J| exp f C{s)dsj (B2) 
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Then, each chronological exponential is replaced by a usual exponential, and the limit n — s> oo 
of the product should converge towards a well defined limit. 

Texp (- [ djC{s)) = lim e""^" ■ e"^-! ■ . . . ■ e"^^ ■ e"^^ (B3) 
\ Jo J 

where Ai is simply 

(piAt/n \ 
- / dsC{s) (B4) 

J{i-l)At/n J 

Thus, the product (IBip can be rewritten as 

lim e^i+^2+-+^"-i+^" ■ e"^" ■ e"^-! ■ . . . ■ e"^^ ■ e"^^ (B5) 

n— >oo 

The last exponential can be further reduced by using the BCH formula, keeping the triple 
commutators: 

exp(A) exp(i?) = exp(^A + B+ ^-[A, B] + ^[A, [A, B]] + ^[B, [5, A]] . . )j (B6) 

Using BCH again on the right hand side of (1B6I) . we find a reverse formula, exact up to 
triple commutators, 

exp(A + B)= exp A] - ^[A, [A, B]] + ^[B, [B, A]]^ exp(A) exp(i?) (B7) 

To reduce the last term in (IBSp . we define the sums 

Si = 0; 

Si = Y,^j = ^i + ^2 + ... + A_,, (B8) 
i=i 

and we apply (1B7P recursively to 



e 2 



g ['S'ti i-Att,] [5*71 ) ['S'n ,-A-n,]] + ^ [Ajt, , [-An ,57;,]] 



= (.-k^7=i[Si,A.]-^j:U[s„[Si,A,]]+lJ:i^,[A„[A,,s,]]+Rn . ■ e^^ ■ . . . e^". (B9) 

In the large n limit, A^ ~ v^At/n, S*, ~ v^At and i?„ ~ 1/n, it results that the term 
I X]r=i[^«' [^«' '^j]] remaining term i?„ do not contribute. Equation fIBSP reduces to 

/ n n \ 

Jim exp -- J2[Si, A] - 3 Et^- ^«]] ' (^^O) 
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or, using Riemann sums instead, 



exp 



At 



ds / du[C{u),C{s)] 



1 



At 



ds I du J dv[C{v),[C{u),Cis)]]] . (Bll) 



When specifying to the Langevin dynamics 

p d 



ds C{s) = ds 



d d 
m or op op 



we obtain 



[ds C{s) , du C{u)] = {tp{u)duds — ■ip{s)dsdu) 



1 d d 
7- 



(B12) 



(B13) 



mdr dp ^ 

which contributes to a At'^/^ term, while the triple commutator gives only a At^/^ term and 



can be disregarded. Finally, 



At 



ds / du[t{u),t{s)] 



At 



s^{s)ds — ds ^p{u) du) 



mdr ^ dp 



'At f^* r 

— / i/?(n)dn— / ds I ip{u)du 
. 2 Jo Jo Jo 



I d _ d 

mdr ^ dp 



(B14) 



and the justification of fl2^ is complete. 

The Magnus expansion cited in eq. (14) of Ref.— is a commutator expansion for the 
logarithm of a chronological exponential: 

(B15) 



n{t) = In (^Texp ( / dsC{s] 
The three first terms of the expansion are 



dsC{s) + - I ds, I ds2[C{s,),C{s^)] 



4 Jo 



S2 



dSl / dS2 / dS3[C{si),[Cis2),Cis3)]] 



Sl 



Sl 



-- / dsi / ds2 I ds,[[C{s,),C{s2)],C{s,)]. 



(B16) 



Using BCH, we get 



exp{fl(t)) = exp ( / ds£(s) ) ■ exp ( - / dsi / ds2[C{si), C{s2)] 



Sl 



S2 



•exp( + - / dsi / ds2 I ds3[C{si),[C{s2),C{s3)]] 



Sl 



Sl 



+ - I ds, I ds2 I dss[[C{s,),C{s2)],C{ss)] 



1 

+ 2 



Sl 



dsi / ds2 [£(si),£(s2)], / dssCiss 



(BIT) 
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Finally, the last commutator reads 



- f dSi r dS2 [ dS3[[CiSi),C{s2)],C{s^)] 
4 -'0 ^0 ^0 







4 Jo Jo Jo 

+ - f ds, r dS2 [ ds^[[C{Si),C{s2)],C{s3)] 
4 In Jo Jsi 



'0 

4 



dSi ['dS2 [ ' dS3[[C{si),jC{s2)]^jC{ss)] 

Jo Jo 

^ ^ ds, Hds, [ ds2[£(s3),[£(si),£(s2)]]. (B18) 

^0 Jsi 



The last term in flBlSp cancels the second line of flB16p . and one gets, after renaming the 
dummy integration variables: 



exp(n(t)) = exp ( I ds C{s) j ■ exp J ds J du [C{s) , C{u)] 



ds du dv[[C{s),Ciu)],Civ)] 
-J Jo Jo Jo 



(B19) 



or equivalently, 

exp(— r2(t)) = exp 



ds£(s)^ ■ exp ^ — - J ds J du [C{u) , C{s)] 



-- ds du dv[C{v),[C{u),C{s)]] 
■J Jo Jo Jo 



(B20) 



identical to flBUj) 



APPENDIX C: THE VELOCITY- VERLET AND LANGEVIN ALGORITHMS 



Let us start with the splitting flTSl) . Each operator acts recursively on the coordinates 
r,p of p{r,p). This leads to the sequence: 

At 



exp{-A/2) p{r,p) = Py-P^^P 
exp{—B)p{r,p) = p {r,p — f{r)At) 
exp{-B) exp{-A/2) p{r,p) = p (r - {p - f{r)At)^,p - f{r)At 



exp(-A/2)exp(-5)exp(-A/2) p{r,p) 

At 



p + p- Atf\ r-p^ j 
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The arguments of p must be identified to the initial position and impulsion {r{t),p{t)) of a 
trajectory ending at the final position r = r{t + At), p = p{t + At). Thus, 



p{t) = pit + At)-f (^r{t + At) - p{t + At)^^ At; 

r{t) = r{t + At) - {p{t + At)+p{t))^. 

2m 



(C2) 

(C3) 



Then, equations can be reversed to express {r{t + At),p(t + At)) as functions of (r(t),p(t)), 
and in O, r(t + At) - p{t + At)At/2m can be replaced by r(t) + p{t)At/2m. 



p{t + At) = p{t) + / (^r(t) + P{t)^^ At; 
r(t + At) = r(t) + (p(t) +p(t + At))^. 

^ III 



Meanwhile, the splitting (|T9i) leads to 



Pit + At) = Pit) + ^ (^/(r(t)) + /(r(t + At)) j 

At At^ 

r(t + At) = rit)+pit)— + firit)) — . 

m 2m 

For the Langevin dynamics of one dimensional systems, we use: 

"At /■*+^* 
— J i^iu)du 

ds / i/'(n)dn 



exp(-A/2) exp(-5) exp(-A/2) ■ exp |cr 



1 d 



1- 



d 



m dr dp 

Swith A = Atp^; B = [At (/(r) - -fp) + a H J- Given that 



exp 



we find 



-A/2 



{Xir) + Yp)^ 



At 



p(r, p) = p ( r,pe^ + X(r) 



p(r,p) = p r-p—,p 
2m 



e ^p(r, p) = p I r, p — I At/(r) + a 



ipiu) du 



-A/2-B 



e '^/^p(r, p) = p ( r — 



At 
2m 



- At/(r~) + a 



p + pe^^* - ( At/(r) + (T V'l^^) dn 

e7Af _ 1 



3^^* - 1 

7At 

■t+At 



(C4) 
(C5) 



(C6) 
(C7) 



p{r,p), (C8) 



(C9) 



37At _ 1 



7At 



t+At 



ipiu) du 



7At 



(CIO) 
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In the above expression, r must be identified with r{t + At), and p with p{t + At). The 



resulting function is e 



~^^'^p{r,p) = p{r',p'), where r' must be replaced with r(t), p' 



with p{t), and r', given by the above expression, with f ^ r(t) +p(t)/2m at the desired 
level of accuracy. These relations can be reverted in such a way that p{t + At) and r(t + At) 
depend explicitely on p(t) and r(t). The result is 

1 - e-^^* 



p{t + At) = p(t)e-^^* + 



At \ /■*+^* 
Atf{r{t)+p{t)—j +a il^{u)du 



7At 



At 



r(t + At) = r(t) + — p{t) + p(t + At) 
2m 



(Cll) 



Finally, the commutator correction simply shifts the previous result by a At'^/^ amount. 



pit + At) = p(t)e-^^* + 
"At 



At/(r(t)+p(t)^^ 



1 - e- 



7At 



+70- 



At 



jp{u)du~ j ds il^{u)du 
Jo Jt 

r{t + At) = r(t) + ( p{t) + p{t + At) ) ^ 



(C12) 



+ - 



m 



At 



ds '4^{u) du 



2m 
At '■*+^* 



■j/>(n) dn 



(C13) 



Alternatively, the other splitting leads to slightly more complex expressions: 



e-''/'e-^e-''/'p{r,p) 



At 

p\r 

m 



pe 



7At/2 



At/(r) + a 



t+At 



■j/?(m) du 



^jAt/2 _ I 



jAt 



pe 



yAt 



rt+At \ ' 

At/(r) + a il){u) du\e^^''^- 



7At/2 _ ]_ 



7At 



t+At \ g7At/2 _ I 

Atfir) + a I duj 



(C14) 



This finally leads to the algorithm 

p{t + At) = p(t)e-^^* + 7-^(1 - e-^^*/2) ( f^^^^ ^ ^ y.^^(^))g-7AV2 



"t+At _ g-7At 

-(7 / ip(u)du- 



7At 



r(t + At) 





'At 


+70- 




r(t)- 


At 


F — 


m 


cr 


- pAt 


H 




m 


.Jo 



-t+At 



p{t)e 



"At 



(C15) 



~-yAt/2 



+ 



1 -e 



--fAt/2 



t+At 



il^{u) du 



r-At PS pt+At 

/ ds ipiu) du — / 'tp{u) du 

Jo Jt ^ Jt 



(C16) 
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These results generalize straightforwardly to the three dimensional many particles case, by 
adding the space and particles indices. For instance 

1 - e-^^* 

It 



At \ /•*+^* 
AtU(r{t)+pit) — j +a I ^p^Mdu 



-fAt 



-70" 



At '•*+^* 



At 



rpiaiu)du- / ds xl^ia{u)du 
't Jo Jt 

Vi^it + At) = ri^{t) + (^Pia{t) + pi^{t + At)^ ^ 

a r z"^* r At f^+^^ 

H — / ds / i^iMdu-— / i^iMdu 
"^VJo Jt ^ Jt 



(C17) 



(C18) 



corresponds to the splitting P3|) . If one expands in powers of At up to At^, the algorithm 
reads 

Pi^it + At) = piait) (^1 - 7At + -^^^ + fia (^r{t) + P{t)^ At (^1 - ^ 



pt+At pt+At PS 

+0- / Via(ii) dn - 0-7 / ds tl^iaiu) du; 
Jt Jt Jt 



(C19) 



It Jt Jt 

/ X At / 7At\ , , , At^ 
ri„(t + At) = ri„(t) +p,,(t)— 1 - ^ + /^a(r(t)) — 



a 

+ — 

m 



t+At 



ds ipia{u) du, 



(C20) 



which are eq (147|) and (l46l) . Alternatively, the splitting (jH]) leads to: 

pUt + At) = pUt)e~^^' + 7-^(1 - e"^^*/2) (^/,,(r(t + At)) + /..(r(t))e-^^*/2^ 

rt+At ^ _ g-7At 

+0- / '4>ic^{u) du- 



7At 



-7a 



At 



t+At 



At 



'il^ia{u)du- I ds il^ic,{u)du 



(C21) 



/ N At 

r,„(t + At) = ri^{t) + — 



+ - 



m 



m 

"At 



Pia{t)e 



-7Ai/2 



+ 



1 -e 



Atfia{r{t)) + a 'i/^ia{u)du 
lAt V 



/At ;-t+At 
ds 'il^ia{u) du - — il)ia{u)du 



(C22) 



and finally, 

/ -i'^At^\ At At^ 

pUt + At) = pUt)\l - 7At + j + — (/,,(r(t)) + fUrit + At)) - 7— /,„(r(t)) 

/t+At pt+At PS 

ipiaiu) dw - 70- / ds / il^ia{u) du] (C23) 

At / At\ At^ cr /■*+^* 

r,,(t + At) = r,,(t) +p.,a—{ 1 - 7^ + fiairit))— + — ds Via(M) dM,(C24) 

m \ 2 J 2m m Jt Jt 



which are eq (|5^ and (l55l) 
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APPENDIX D: DERIVATION OF THE DPD ALGORITHM 



The operators exp{—A) and exp{—B), with A and B introduced in (!66|) . act respectively 
on test functions hke 



e-^p(r,p) = p(r,e^*^(") ■p-r(r)^^(e^*^(") - 1) ■K(r;t) 

e~^p{r,p) = p{r p,p). 

m 



(Dl) 
(D2) 



The sphtting e "^^^e ^1"^ leads to, following manipulations similar to appendix ( ICl) : 

At 



r{t + At) = r{t) + l^p(t) + p(t + At) 
p(t + At) = e-^*^(^W+PWi^)p(t) 



2m' 



(D3) 



+r(r(t)+p(t)^) (l-e-^*^W*HPW^))-K(r(t)+p(t)^;t), (D4) 

which, upon expansion in powers of At up to At^, gives the algorithm (ITTil and (ITOll . The 
second splitting e~^/^e~'^e~'^/^ leads to: 



r(t + At) = r(t) H 

m 



;(D5) 



p(t + At) 



g-^r(r{t+Ai))g-^r{r(t)) 



e-^r(r.(t))p^^) ^ r(r(t))-^ (^1 - e-^^^^*)) j • K(r(t); t) 
+e-^^W*+^*«r(r(t))-i(^l - e-^^W*))^ ■ K(r(t);t) 

+ r(r(t + At))-^ (^1 - r(r(t+At))^ . ^(^(^ _^ 

which, upon expansion in powers of At up to At^, gives the algorithm (!73l) and fl72|) . The 
correction term [i?, [A, i?]] acts only on the impulsions: 



\B\A,B\\ = 2- At V 



ijkl,af3 



t+At 



X 



4jij{u)du 
d (w{rki) 



t+At 



{rkf3 - rip) 



%l)ki{u)<lu 
d 



(D7) 



As a further simplification, we remark that the terms involving two different pairs ij and kl 



are of zero mean, uncorrelated from any 



J^~^^^ipjrin{u)du , and of variance At^. Therefore, 



within a weak order two scheme, we can replace (]D7[) by its average value and the correction 
term becomes: 



[B,[A,B]] = 2-At'Y, 



ij,af3 



d 



d 



, (D8) 
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and finally, each impulsion is incremented by an amount 



w{ri,j) 



mc 



d (w{rij) 



(D9) 



where c is 12 for the splitting e ^/"^e ^1"^ or -24 for the splitting e ^I'^e "^Z^. 

The correction exp(— 1/2 JJdnds [£(u),£(s)]) gives terms proportional to 



J(s — At/2)^/;jj(s)ds, which have been shown to be uncorrelated to the J^^^^ipij{u)du 
vanishing mean value and of variance At^. These terms are the pendant of the corrections 
in the Ermak algorithm, and can be neglected in practice. 

The correction term exp(— 1/3 jjj dvduds [jC{v), [jC{u) , jC{s)]]) can be written as 



of 



cr 

3m 



i,j,k,l 



tlJij{u)du ipki{s){s - t)ds + / ipki{u)du 'iljij{s){s - t)ds 



- 2ds^j^ ^pij{u)du j ^ki{u)d 



a,/3 



or, 



rki 



dpkp j ' 



If the pairs and {k,l) differ, the random increment turns out to be uncorrelated to 



the 



u 



of variance At^ and therefore can be evacuated from an order two 
integration scheme. However, the term corresponding to the same pairs gives an positive 
At^ contribution (we use here the Stratonovitch interpretation of the stochastic integrals). 



t+At 



'- Jt 



ipij{u)du iljij{s){s — t)ds — ds J ipij{u)du 
— 1^ / %j[u)du - - I ds 



4jij{u)du 



-At' 



(DIO) 



while the last term of equation (IDlOp reduces to (IDSp . leading to an operator: 



3m 



t+At 



a,l3 



i/jij{u)du ipij{s){s — t)ds — ds / ilJij{u)du 



1 2 



w{rij) 



d /w(r. 



l^ia ^ ja) ^ 



13) 



d 



d 



dpi(3 dpjji 



(Dll) 



and the corresponding corrective term for the impulsions Pis equals 



{Ap 



3m 



w{rij) 



(D12) 
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Together, the correction of the triple commutators is the sum of {^pip)i (1D9P and 
{^Pi0)2 are 



— e - 3 |: Ur'^- - ^-'*-(^(^-^ - ^'"O ■ (^^^' 

For a fixed couple of indices i and j, the differentiation takes the form X ■ VX, with X = 
[w(rjj)/rjj]rjj. A further simplification comes, using the identity 
X ■ VX = V(X2)/2 - X A (V A X), 

V. A (^(.. - = V. {^^) A (n,) + (^) A (V. A n,) 

= 0. (D14) 

Finally, the correction reduces to 

Surprisingly, the term cancels out if c = 12, meaning that there is no contribution at all in 
the case of the splitting e~'^/^e~^e~"^/^, and that the algorithm (!70| [7T|) is already a genuine 
weak order two algorithm for DPD. 

However, the correction term must be accounted for when the splitting e~^l'^e~^e~^l'^ is 
considered, ending with the weak order two algorithm (!77 l [78l 1791) . 
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CAPTIONS 



FIG 1. Illustration of the Fokker-Planck evolution in phase space. 

FIG 2. Respective efficiencies of the Ermak and Ricci-Ciccotti algorithms. See text for 
details. 

FIG 3. Pair correlation function g{r) for the ideal gas and a standard weight function 
w{r) = 1 — r/vc- From left to right, At takes the values 0.01, 0.05 and 0.1. The algorithms 
presented here are ABA, BAB, BAB deprived of its correcting term, Shardlow SI and 
DPD-VV. Simulations involve 4000 particles, a run length of 60000 time steps, and 
parameters values 7 = 4.5, T = 1., a = 3. 

FIG 4. Thermalization properties of ABA, BAB, Shardlow SI, DPD-VV (singular 
weight w{r)), ABA and DPD-VV (regular weight wi{r)). On the left panel: simulated 
temperatures (T) vs time step At; on the right panel: a semi-log plot of the difference 
I (T) — 1\ vs the time step At. 

FIG 5. Pair correlation function g{r) for an interacting fluid, a weight function 
w{r) = 1 — r/r^ and a pair interaction potential aw^(r)/2, with a = 25. From left to 
right, we show the algorithm Shardlow SI, ABA and DPD-VV. The full line represents the 
result for At = 0.01, the dashed line for At = 0.05, and the symbols o for At = 0.1. All 
three algorithms coincide in the limit At = 0.01. The ABA algorithm is more accurate for 
At = 0.1 than the two others, a feature that might be explained by its order two nature. 
In the presence of repulsive interaction, the singularity of w{r) has no consequences. 

FIG 6. Presentation of the DPD diffusion coefficient D{At) for an interacting fluid, for 
three different values of the time step At. Left panel: For each algorithm, the curves are 
normalized by their value at At = 0, 01, and only relative variations of D{At) are shown, as 
in Ref.-. Right panel: the bare D{At) are shown. Depending on the algorithm chosen, the 
self- diffusion coefficient changes by a proportion of 6%, which is also our estimate of the 
statistical error affecting our data. We conclude that the values obtained with At = 0.05 
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are consistent with the reference values obtained with At = 0.01, and that there is no 
significant difference between these three algorithms, regarding the self-diffusion properties. 
For completeness, we calculated the self-diffusion coefficient for the ideal gas, and found 
that the data obtained for the three algorithms were almost indistinguishable (all within 
0.5%). The self-diffusion coefficient of the ideal gas is close to 0.535. The interaction 
potential reduces only lightly this value, excluding any caging effect or thermal activation. 

FIG 7. Pair correlation function g{r) for the ideal gas and a weight function 
Wi{r) = (r/rc) x (1 — r/r^. From left to right. At takes the values 0.01, 0.05 and 0.1. 
The algorithms presented here are ABA, Shardlow SI and DPD-VV. Simulations involve 
4000 particles, a run length of 60000 time steps, and parameters values 7 = 4.5, T = 1., 
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